##########################
#Alternative Explanations#
##########################
#This to obtain:
#table_A5.tex
#table_A6.tex

library("broom")
library("dplyr")
library("grid")
library("tidyr")
library("sandwich")
library("lfe")
library("broom")
library("dplyr")
library("tidyr")
library("grid")
library("stargazer")
library("lmtest")
library("multiwayvcov")
library("ggplot2")
library("glue")
library("spdep")
library("spatialreg")
library("stringr")
library("readxl")  
library("rgdal")

setwd("/Users/bgpopescu/")
#setwd("C:/Users/bogdanp/")

#############
#1921 Census#
#############

census_1921 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                       layer="census_1921_points")


census_1921$pct_serbs_1921<-census_1921$Christian_orthodox/census_1921$Total*100
census_1921$pct_croats_1921<-census_1921$Roman_catholic/census_1921$Total*100

census_1921$dist1 <- as.numeric( with (census_1921,ifelse(census_1921$treat==1, 1, -1)))
census_1921$dist2 <-census_1921$dist1*(census_1921$krajna6_distance/1000)

census_1921$quad1<-ifelse((census_1921$lon > 18 & census_1921$lon< 20) & 
                            (census_1921$lat<47  & census_1921$lat>44), 1, 0)

census_1921$quad2<-ifelse((census_1921$lon > 16 & census_1921$lon< 18) & 
                            (census_1921$lat<47  & census_1921$lat>44), 1, 0)

census_1921$quad3<-ifelse((census_1921$lon > 14 & census_1921$lon< 16) & 
                            (census_1921$lat<47  & census_1921$lat>44), 1, 0)

census_1921$bfe1 <- ifelse(census_1921$krajna6_NEAR_FID == 1, 1,0)
census_1921$bfe2 <- ifelse(census_1921$krajna6_NEAR_FID == 2, 1,0)
census_1921$bfe3 <- ifelse(census_1921$krajna6_NEAR_FID == 3, 1,0)
census_1921$bfe4 <- ifelse(census_1921$krajna6_NEAR_FID == 4, 1,0)
census_1921$bfe5 <- ifelse(census_1921$krajna6_NEAR_FID == 5, 1,0)
census_1921$bfe6 <- ifelse(census_1921$krajna6_NEAR_FID == 6, 1,0)
census_1921$bfe7 <- ifelse(census_1921$krajna6_NEAR_FID == 7, 1,0)
census_1921$bfe8 <- ifelse(census_1921$krajna6_NEAR_FID == 8, 1,0)


census_1921$POINT_X<-census_1921$lon
census_1921$POINT_Y<-census_1921$lat
census_1921<-subset(census_1921, !is.na(treat))


#############
#1931 Census#
#############

census_1931 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                       layer="census_1931_points")

census_1931$pct_serbs_1931<-census_1931$Orthodox/census_1931$total*100
census_1931$pct_croats_1931<-census_1931$Roman_catholic/census_1931$total*100
census_1931$dist1 <- as.numeric( with (census_1931,ifelse(census_1931$treat==1, 1, -1)))
census_1931$dist2 <-census_1931$dist1*(census_1931$krajna6_distance/1000)
census_1931$krajna6_NEAR_FID

census_1931$quad1<-ifelse((census_1931$lon > 18 & census_1931$lon< 20) & 
                      (census_1931$lat<47  & census_1931$lat>44), 1, 0)

census_1931$quad2<-ifelse((census_1931$lon > 16 & census_1931$lon< 18) & 
                      (census_1931$lat<47  & census_1931$lat>44), 1, 0)

census_1931$quad3<-ifelse((census_1931$lon > 14 & census_1931$lon< 16) & 
                      (census_1931$lat<47  & census_1931$lat>44), 1, 0)


census_1931$bfe1 <- ifelse(census_1931$krajna6_NEAR_FID == 1, 1,0)
census_1931$bfe2 <- ifelse(census_1931$krajna6_NEAR_FID == 2, 1,0)
census_1931$bfe3 <- ifelse(census_1931$krajna6_NEAR_FID == 3, 1,0)
census_1931$bfe4 <- ifelse(census_1931$krajna6_NEAR_FID == 4, 1,0)
census_1931$bfe5 <- ifelse(census_1931$krajna6_NEAR_FID == 5, 1,0)
census_1931$bfe6 <- ifelse(census_1931$krajna6_NEAR_FID == 6, 1,0)
census_1931$bfe7 <- ifelse(census_1931$krajna6_NEAR_FID == 7, 1,0)
census_1931$bfe8 <- ifelse(census_1931$krajna6_NEAR_FID == 8, 1,0)

census_1931$POINT_X<-census_1931$lon
census_1931$POINT_Y<-census_1931$lat

census_1931<-subset(census_1931, !is.na(treat))


###########################
#1991, 2001, 2011 Censuses#
###########################
data <- read_excel("./Dropbox/Legacies_Project/Analysis/data/merge.xlsx", sheet=1, col_names = TRUE, skip = 0)
data$treat <- as.numeric( with (data,ifelse((data$treat_KRAJ_kraj6==1), 1, 0)))
data$dist1 <- as.numeric( with (data,ifelse(data$treat==1, 1, -1)))
data$dist2 <-data$dist1*data$krajna6_distance


HRV_adm3 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm3")
HRV_adm3<-st_simplify(HRV_adm3, dTolerance = 0.005)
HRV_adm3<-subset(HRV_adm3, select = c(ID_2))
HRV_adm3_data = read_excel("./Dropbox/Legacies_Project/Analysis/data/merge.xlsx", sheet=1, col_names = TRUE, skip = 0)
data<-left_join(HRV_adm3, HRV_adm3_data, by = c("ID_2"="ID_2"))


##########
#Analysis#
##########

pct_serbs_1921<-lm(pct_serbs_1921~treat + POINT_X + POINT_Y + zagreb_distance +
         bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
         quad1 + quad2 + quad3, data = census_1921)


pct_croats_1921<-lm(pct_croats_1921~treat + POINT_X + POINT_Y + zagreb_distance +
         bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
         quad1 + quad2 + quad3, data = census_1921)


pct_serbs_1931<-lm(pct_serbs_1931~treat + POINT_X + POINT_Y + zagreb_distance +
         bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
         quad1 + quad2 + quad3, data = census_1931)

pct_croats_1931<-lm(pct_croats_1931~treat + POINT_X + POINT_Y + zagreb_distance +
                     bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
                     quad1 + quad2 + quad3, data = census_1931)


pct_serbs_1931<-lm(pct_serbs_1931~treat + POINT_X + POINT_Y + zagreb_distance +
                      bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
                      quad1 + quad2 + quad3, data = census_1931)

pct_serbs_1991<-lm(pct_serbs_1991~treat + POINT_X + POINT_Y + zagreb_distance +
         bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
         quad1 + quad2 + quad3, data = data)


pct_croats_1991<-lm(pct_croats_1991~treat + POINT_X + POINT_Y + zagreb_distance +
                     bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
                     quad1 + quad2 + quad3, data = data)


pct_serbs_1991<-lm(pct_serbs_1991~treat + POINT_X + POINT_Y + zagreb_distance +
                      bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
                      quad1 + quad2 + quad3, data = data)


pct_croats_2001<-lm(Pct_Croats_2001~treat + POINT_X + POINT_Y + zagreb_distance +
                      bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
                      quad1 + quad2 + quad3, data = data)


pct_serbs_2001<-lm(Pct_Serbs_2001~treat + POINT_X + POINT_Y + zagreb_distance +
                     bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
                     quad1 + quad2 + quad3, data = data)


pct_change_croats<-lm(pct_change_croats~treat + POINT_X + POINT_Y + zagreb_distance +
                     bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
                     quad1 + quad2 + quad3, data = data)

pct_change_serbs<-lm(pct_change_croats~treat + POINT_X + POINT_Y + zagreb_distance +
                        bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
                        quad1 + quad2 + quad3, data = data)



pct_croats_2011<-lm(Pct_Croats~treat + POINT_X + POINT_Y + zagreb_distance +
                      bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
                      quad1 + quad2 + quad3, data = data)


pct_serbs_2011<-lm(Pct_Serbs~treat + POINT_X + POINT_Y + zagreb_distance +
                     bfe1 + bfe2 + bfe3 + bfe4 + bfe5 + bfe6 + bfe7 + 
                     quad1 + quad2 + quad3, data = data)



######################
#Step1: RELEVANT VARS#
######################

relevant_vars<-c("pct_serbs_1921",
                 "pct_serbs_1931",
                 "pct_serbs_1991",
                 "Pct_Serbs_2001",
                 "pct_change_serbs",
                 "Pct_Serbs")


#####################
#Step2: MEANS AND SD#
#####################

mean_pct_serbs_1921<-round(mean(census_1921$pct_serbs_1921, na.rm=T), digits = 3)
mean_pct_serbs_1931<-round(mean(census_1931$pct_serbs_1931, na.rm=T), digits = 3)
mean_pct_serbs_1991<-round(mean(data$pct_serbs_1991, na.rm=T), digits = 3)
mean_Pct_Serbs_2001<-round(mean(data$Pct_Serbs_2001, na.rm=T), digits = 3)
mean_pct_change_serbs<-round(mean(data$pct_change_serbs, na.rm=T), digits = 3)
mean_Pct_Serbs<-round(mean(data$Pct_Serbs, na.rm=T), digits = 3)
means<-c(mean_pct_serbs_1921,
         mean_pct_serbs_1931,
         mean_pct_serbs_1991,
         mean_Pct_Serbs_2001,
         mean_pct_change_serbs,
         mean_Pct_Serbs)


sd_pct_serbs_1921<-round(sd(census_1921$pct_serbs_1921, na.rm=T), digits = 3)
sd_pct_serbs_1931<-round(sd(census_1931$pct_serbs_1931, na.rm=T), digits = 3)
sd_pct_serbs_1991<-round(sd(data$pct_serbs_1991, na.rm=T), digits = 3)
sd_Pct_Serbs_2001<-round(sd(data$Pct_Serbs_2001, na.rm=T), digits = 3)
sd_pct_change_serbs<-round(sd(data$pct_change_serbs, na.rm=T), digits = 3)
sd_Pct_Serbs<-round(sd(data$Pct_Serbs, na.rm=T), digits = 3)
sds<-c(sd_pct_serbs_1921,
       sd_pct_serbs_1931,
       sd_pct_serbs_1991,
       sd_Pct_Serbs_2001,
       sd_pct_change_serbs,
       sd_Pct_Serbs)



####################
#Step3: LIST MODELS#
####################

list_models<-list(pct_serbs_1921,
                  pct_serbs_1931,
                  pct_serbs_1991,
                  pct_serbs_2001,
                  pct_change_serbs,
                  pct_serbs_2011)


######################
#Step5: Column labels#
######################


column_labels= c("1921",
                 "1931",
                 "1991",
                 "2001",
                 "Change '91-'01",
                 "2011")



serbs<-stargazer(list_models,
                        column.labels= column_labels,
                        keep = c("treat"),
                se = list(coeftest(pct_serbs_1921, cluster.vcov(pct_serbs_1921, census_1921$Kotar_Municipality))[, 2],
                          coeftest(pct_serbs_1931, cluster.vcov(pct_serbs_1931, census_1931$srez_name))[, 2],
                          coeftest(pct_serbs_1991, cluster.vcov(pct_serbs_1991, data$NAME_2))[, 2],
                          coeftest(pct_serbs_2001, cluster.vcov(pct_serbs_1991, data$NAME_2))[, 2],
                          coeftest(pct_change_croats, cluster.vcov(pct_change_croats, data$NAME_2))[, 2],
                          coeftest(pct_croats_2011, cluster.vcov(pct_croats_2011, data$NAME_2))[, 2]),
                covariate.labels=c("Military Territory"),
                dep.var.caption = "Dependent variable",
                dep.var.labels.include = F,
                omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                add.lines=list(c("Mean", means),
                               c("SD", sds),
                               c("Boundary FE", rep("Yes", 7)),
                               c("Region FE", rep("Yes", 7))))
note_latex <- "\\multicolumn{7}{l} {\\parbox[t]{16cm}{ \\footnotesize \\textit{Notes:} 
Coefficients and robust standard errors in parantheses from OLS regression. 
$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.
The units of analysis are settlements for the corresponding censuses.
See codebook in the appendix for data sources and how the variables were calculated.}}"
serbs[grepl("Note",serbs)] <- note_latex
cat(serbs, file="./Dropbox/Legacies_Project/Paper/tables/table_A5.tex", sep="\n")


######################
#Step1: RELEVANT VARS#
######################


relevant_vars<-c("pct_croats_1921",
                 "pct_croats_1931",
                 "pct_croats_1991",
                 "Pct_Croats_2001",
                 "pct_change_croats",
                 "Pct_croats")

#####################
#Step2: MEANS AND SD#
#####################

mean_pct_croats_1921<-round(mean(census_1921$pct_croats_1921, na.rm=T), digits = 3)
mean_pct_croats_1931<-round(mean(census_1931$pct_croats_1931, na.rm=T), digits = 3)
mean_pct_croats_1991<-round(mean(data$pct_croats_1991, na.rm=T), digits = 3)
mean_pct_croats_2001<-round(mean(data$Pct_Croats_2001, na.rm=T), digits = 3)
mean_pct_change_croats<-round(mean(data$pct_change_croats, na.rm=T), digits = 3)
mean_pct_croats<-round(mean(data$Pct_Croats, na.rm=T), digits = 3)
means<-c(mean_pct_croats_1921,
         mean_pct_croats_1931,
         mean_pct_croats_1991,
         mean_pct_croats_2001,
         mean_pct_change_croats,
         mean_pct_croats)


sd_pct_croats_1921<-round(sd(census_1921$pct_croats_1921, na.rm=T), digits = 3)
sd_pct_croats_1931<-round(sd(census_1931$pct_croats_1931, na.rm=T), digits = 3)
sd_pct_croats_1991<-round(sd(data$pct_croats_1991, na.rm=T), digits = 3)
sd_pct_croats_2001<-round(sd(data$Pct_Croats_2001, na.rm=T), digits = 3)
sd_pct_change_croats<-round(sd(data$pct_change_croats, na.rm=T), digits = 3)
sd_pct_croats<-round(sd(data$Pct_Croats, na.rm=T), digits = 3)
sds<-c(sd_pct_croats_1921,
       sd_pct_croats_1931,
       sd_pct_croats_1991,
       sd_pct_croats_2001,
       sd_pct_change_croats,
       sd_pct_croats)


####################
#Step3: LIST MODELS#
####################

list_models<-list(pct_croats_1921,
                  pct_croats_1931,
                  pct_croats_1991,
                  pct_croats_2001,
                  pct_change_croats,
                  pct_croats_2011)


croats<-stargazer(list_models,
                  column.labels= column_labels,
                  keep = c("treat"),
                  se = list(coeftest(pct_croats_1921, cluster.vcov(pct_croats_1921, census_1921$Kotar_Municipality))[, 2],
                            coeftest(pct_croats_1931, cluster.vcov(pct_croats_1931, census_1931$srez_name))[, 2],
                            coeftest(pct_croats_1991, cluster.vcov(pct_croats_1991, data$NAME_2))[, 2],
                            coeftest(pct_croats_2001, cluster.vcov(pct_croats_1991, data$NAME_2))[, 2],
                            coeftest(pct_change_croats, cluster.vcov(pct_change_croats, data$NAME_2))[, 2],
                            coeftest(pct_croats_2011, cluster.vcov(pct_croats_2011, data$NAME_2))[, 2]),
                  covariate.labels=c("Military Territory"),
                  dep.var.caption = "Dependent variable",
                  dep.var.labels.include = F,
                  omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                  add.lines=list(c("Mean", means),
                                 c("SD", sds),
                                 c("Boundary FE", rep("Yes", 7)),
                                 c("Region FE", rep("Yes", 7))))
note_latex <- "\\multicolumn{7}{l} {\\parbox[t]{16cm}{ \\footnotesize \\textit{Notes:} 
Coefficients and robust standard errors in parantheses from OLS regression. 
$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01.
The units of analysis are settlements for the corresponding censuses. 
See codebook in the appendix for data sources and how the variables were calculated.}}"
croats[grepl("Note", croats)] <- note_latex
cat(croats, file="./Dropbox/Legacies_Project/Paper/tables/table_A6.tex", sep="\n")

